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A mechanism for fast magnetic reconnection in collisionless plasma is studied for understand- 
ing sawtooth collapse in tokamak discharges. Nonlinear growth of the tearing mode driven 
by electron inertia is analytically estimated by invoking the energy principle for the first time. 
Decrease of potential energy in the nonlinear regime (where the island width exceeds the elec- 
tron skin depth) is found to be steeper than in the linear regime, resulting in acceleration of the 
reconnection. Release of free energy by such ideal fluid motion leads to unsteady and strong 
convective flow, which theoretically corroborates the inertia-driven collapse model of the saw- 
tooth crash [D. Biskamp and J. F. Drake, Phys. Rev. Lett. 73, 971 (1994)]. 

1 Introduction 

Sawtooth collapse in tokamak plasmas has been a puzzling phenomena for decades. Although 
the m = 1 kink-tearing mode is essential for onset of this dynamics, Kadomtsev's full reconnec- 
tion model HI and nonlinear growth of the resistive m = 1 mode [2] (both based on resistive 
magnetohydrodynamic theory) fails to explain the short collapse times (~ 100/zs) as well as 
partial reconnections observed in experiments. Since resistivity is small in high-temperature 
tokamaks, two-fluid effects are expected to play an important role for triggering fast (or explo- 
sive) magnetic reconnection as in solar flares and magnetospheric substorms. 

In earlier works @|, the linear growth rate of the kink-tearing mode in the collisionless 
regime has been analyzed extensively by using asymptotic matching, which shows an enhance- 
ment of the growth rate due to two-fluid effects, even in the absence of resistivity. Furthermore, 
direct numerical simulations E [6l 13 of two-fluid models show acceleration of reconnection 
in the nonlinear phase, which indicates explosive tendencies until numerical error or artificial 
dissipation terminates them. 

However, theoretical understanding of such explosive phenomena is not yet established due 
to the lack of analytical development. In contrast to the quasi-equilibrium analysis developed 
for resistive reconnections El [8]] , the explosive process of collisionless reconnection should be 
a nonequilibrium problem, in which inertia is not negligible in the force balance and hence 
leads to acceleration of flow. The convenient assumption of steady reconnection is no longer 
appropriate. 

Recent theories (9l [lOl dH emphasize the Hamiltonian nature of two-fluid models and try to 
gain deeper understanding of collisionless reconnection in the ideal limit. 

The purpose of the present work is to predict explosive growth of the kink-tearing mode 
analytically by developing a new approach that is based on the energy principle [fT2l . For sim- 
plicity, we will consider only the effect of electron inertia, which is an attractive mechanism for 
triggering fast reconnection in tokamaks; estimates of the reconnection rate are favorable [fT3lh 
nonlinear acceleration is possible 0, and even the mysterious partial reconnection may be ex- 
plained by an inertia-driven collapse model 031 Q3]]. While we address the same problem as 



Ref. 0, the estimated nonlinear growth is quantitatively different from that of Ref. [6j. Our 
result is confirmed by a direct numerical simulation and its implications for sawtooth collapse 
are discussed in the final section. 



2 Free energy source of tearing induced by electron inertia 



We analyze the following vorticity equation and (collisionless) Ohm's law for velocity field 

v = e z x Vcj)(x, y, t) and magnetic field B = Vip(x, y, t) x e z + B e z : 
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where [/, g] = (V/ x Vg) • e z . The parameter d e denotes the electron skin depth, which is much 
smaller than the system size (d e <C L x ). Since the frozen-in flux for Eq. © is not magnetic 
flux i/; but the electron canonical momentum defined by i(j e = t/j — dlV 2 ^, the effect of electron 
inertia permits magnetic reconnection within a thin layer (~ d e ) despite a lack of resistivity. In 
the same manner as Ref. [6], we consider a static equilibrium state, 



0, ip^ (x) = ipo cos ax, 



(3) 



on a doubly-periodic domain D = [—L x /2, L x /2] x [—L y /2, L y /2] (where a = 2ir/L x ), and 
analyze nonlinear evolution of the tearing mode whose wavenumber in the ^-direction is k = 
2ir/ 'L y at its early linear stage. For sufficiently small k such that 
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this instability is similar to the m = 1 kink-tearing mode in tokamaks (which is marginally 
stable in the ideal MHD limit, d e = 0). FIG. Q] shows contours of ip calculated by direct 
numerical simulation, where e denotes maximum displacement in the x-direction. As shown in 
FIG. [2l the growth of e accelerates when e = e/d e > 1 which is faster than exponential 
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Figure 1: Contours of i/j when e = 4.2rf e 
(d e /L x = 0.01 and L y /L x = 47r) 
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Figure 2: Growth of e = e/d e 
(d e /L x = 0.01 and L y /L x = 47r) 



In order to assess the free energy available from the equilibrium state, we solve the con- 
servation law © for i\) e = tj) — dgV 2 ^ by introducing an incompressible flow map G t : 
D —> D, which depends on time and corresponds to the identity map (G_oo = Id) when 



t = — oo. Let (x,y)(t) = G t (x ,yo) be orbits of fluid elements labeled by their position 
0,2/0) at * = ~~ 00 • Then, the velocity field (or cj)) is related to G t by dG t /dt(x ,yo) = 
e z x V0(x, y, t). Provided that we regard G t as an unstable fluid motion emanating from the 
equilibrium state ©, we can solve © by ij) e {x,y,t) = ip e (Gt(x ,yo) ) t) = ipi°\x ), where 

cos(cke) ~ ij)^\x). By adapting Newcomb's Lagrangian theory lfT6ll . 
we define the Lagrangian for the fluid motion G t as 

L[G t ]=K[G t }-W[G t ], (5) 

where 

K[G t ] = \f \Vct)\ 2 d 2 x and W[G t ] = \ [ (| V^| 2 + d 2 | V 2 ^| 2 ) d 2 x. (6) 

One can confirm that the variational principle 5 J L[G t ]dt = with respect to 5G t yields the 
vorticity equation (OD). 

Note that W plays the role of potential energy and the equilibrium state @ initially stores it 
as free energy. In the same spirit as the energy principle lfT2l . if the potential energy decreases 
(5W < 0) for some displacement map G t , then such a perturbation will grow with the release 
of free energy. 



3 Energy principle for linear stability analysis 



In our linear stability analysis, the equilibrium state is perturbed by an infinitesimal displace- 
ment, G t (xo, yo) = (x , yo) + £( x o> 2/0 3 1)> where £ is a divergence-free vector field on D. We 
seek a linearly unstable tearing mode in the form 
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with a growth rate e(t) oc e 7t . We normalize the eigenfunction £(x) by max \£(x)\ = 1 so 
that e(t) is equal to the maximum displacement in the x-direction and, hence, measures the half 
width of the magnetic island. 

Upon omitting " from equilibrium quantities, ip^ , , , etc., to simplify the notation, 
the eigenvalue problem can be written in the form 

- ~{i 2 /k 2 + €) £'] ' + k 2 ( 7 2 A 2 + €) i = W'i + ^ 2 v 2 T -^v 2 (t//|), (8) 

where V 2 should be interpreted as V 2 = d 2 — k 2 and the prime (0 denotes the x derivative. 
Note, ([8]) ranks as a fourth order ordinary differential equation (unless d e = 0) because of the 
integral operator (1 — <i 2 V 2 ) -1 on the right hand side. By multiplying the both sides of ([8]) by 
£ and integrating over the domain, we get —^ 2 I^ = where 



j(2) = 

w (2) = 



L x /2 



L x /2 
L x /2 

■L x /2 



dx 



w\ 2 - 
- m) 



V 2 



1 - d? e W 2 



(9) 
(10) 



The functionals and are, respectively, related to the kinetic and potential energies 

for the linear perturbation. Hence, by invoking the energy principle |fl"2ll (or the Rayleigh-Ritz 



method), we can search for the most unstable eigenvalue (7 > 0) by minimizing /I^ with 
respect to £. 

Since we assume the ordering © that corresponds to the kink-tearing mode, the eigenfunc- 
tion £ is approximately constant except for thin boundary layers at x = 0, ±L x /2 and has dis- 
continuities around them because of the singular property of ([8]) in the limit of (7/fc) , fc, d e —> 0. 
The electron inertia effect would smooth out these discontinuities. 
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Figure 3: Test function that mimics the un- 
stable tearing mode Fi § ure 4: The linear § rowth rate 7 calcu- 
lated by simulation (d e /L x = 0.01) 

Let us a priori choose the piecewise-linear test function shown in FIG. [H By substituting 
this function into © and ([TOlh we can make negative and keep 1^ finite as follows: 

1^ ~ 4/d e fc 2 , W {2) ~ -2 (1/3 + 9e~ 2 ) d e ^ 2 , where B' y0 = a 2 ip and we have extracted only 
the leading-order term. The linear growth rate is therefore estimated as 

7 = =^/o.776r - 2 = 0.881r -\ (11) 

where r _1 = d e kB' y0 . This result agrees with the general dispersion relation derived by asymp- 
totic matching EH]]. Of course, our analytical estimate of the growth rate depends on how good 
the chosen test function mimics the genuine eigenfunction. Nevertheless, the result predicted 
by the simple function in FIG. [3] shows a satisfactory agreement with the numerically calculated 
growth rate (see FIG. |4]) in the small k region corresponding to the ordering ©. 

4 Variational estimate of explosive nonlinear growth 

Next, we consider the nonlinear phase of the linear instability discussed above. We remark 
in advance that a higher-order perturbation analysis of the Lagrangian (i.e., weakly nonlinear 
analysis) ifTTll will not be successful. Such a perturbation expansion will fail to converge when 
the displacement e (or the island width) reaches the boundary layer width (~ d e ), since the 
eigenfunction has a steep gradient £' ~ ^/d e inside the boundary layers (see FIG. [3]). The 
naive perturbation analysis is, therefore, only valid for < e <C d e9 while e actually exceeds d e 
without saturation as in FIG. El 

To avoid difficulties of a rigorous fully-nonlinear analysis, we again take advantage of the 
variational approach. Namely, we devise a trial fluid motion (parameterized by the amplitude 
e) that tends to decrease the potential energy W as much as possible. When such a motion is 
substituted into the Lagrangian ©, it is expected to be nonlinearly unstable. 

Owing to the symmetry of the mode pattern, it is enough to discuss the boundary layer at 
x = and, moreover, focus on only the 1st quadrant, < x and < y < L y /2. In a heuristic 
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Figure 5: Deformation of contours of ip e by the displacement map (fT2|) 



2xo-^(x ), ^ + f<|/o<^. (iii) 



way based on the simulation result, we consider a displacement map G e : (x , y ) \-> (x, y) 
where the displacement in the x direction is prescribed by 

9e(x ), < y < x - 2' W 

+ f (yo - ^) bo - &(z )], ^ - | < Vo < ^ + i (ii) (12) 

% + I < I/O < ^ 

The regions (i)-(iii) are indicated in FIG. Oleft) and we furthermore define g e as 

{e _e x , < xq < d e , 

d e e^V _1 , d e <x <d e + e, (13) 
x -e, d e + e< x . 

As illustrated in FIG. fright), this displacement map deforms the contours of ip e into a Y-shape. 
From this deformation we find that the potential energy decreases as follows: 



SW[G e ] = - L y B»4 



| + 0(?) 



(14) 



in a nonlinear regime d e <C e <C L x . To obtain the estimate (fT4l) that is likely close to the 
steepest descent, we have technically chosen the map (IT2l) based on the following observations: 

• Around the X points, the flux ip e of the red regions of FIG. Hleft) is squeezed into the 
boundary layers in FIG. fright). On the other hand, the flux is expanded around the O 
points and the areas of the blue regions of FIG.^left) are almost doubled in FIG. fright). 
Since ip e ~ ip except for the boundary layers, both deformations are found to decrease 
magnetic energy (1/2) J \\7ip\ 2 d 2 x as e 3 when d e « e « L^. 

• As is also shown in Ref. 0, a strong current spike develops inside the boundary layers 
[i.e., the red regions in FIG. fright)] in the form of J ~ €B' y0 log \x/d e \ for e = e/d e 3> 1, 
which increases the current energy (1/2) J d 2 J 2 d 2 x (where J = — V 2 ^). However, this 
logarithmic singularity is square-integrable and the current energy change is, at most, of 
the second order 0(e 2 ) in (fT4l) . 



• Only in the intermediate region (ii) located between the X and O points, does the potential 
energy tend to increase. But, we can omit the detailed analysis of this region by taking its 
width I to be sufficiently small: I <C L y . We are allowed to use this approximation as far 
as the kink-tearing ordering © is concerned, in which L y is the longest length scale. 

By introducing time-dependence in e(t), we also need to calculate the kinetic energy, which 
eventually results in 
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where i — t/r . Therefore the Lagrangian © reduces to 
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where U(e) = -(3/tt 2 log2)e 3 + 0(e 2 ) = -0.439e 3 + 0(e 2 ). In the linear regime (e < 1), 
we have already shown that the potential energy decreases as U{e) = — 0.776e 2 . The steeper 
descent where U (e) = — 0.439e 3 in the nonlinear regime (e 3> 1) indicates an explosive growth 
of e during a finite time ^ r . 

We remark that the nonlinear force F(e) = —U'(e) ~ 0{e 2 ) obtained here is different from 
F{e) ~ 0(e 4 ) in the earlier work [6J. While the similar fluid motion around the X and O points 
is considered in Ref. J61, they directly integrate the vorticity equation (OD) over the quadrant 
[0, L x /2] x [0, Ly/2] and arrive at an equation of motion d 2 e/di 2 = F(e) ~ 0(e 4 ). However, 
unless the assumed trial motion happens to be an exact solution, their treatment may lead to a 
wrong equation of motion that does not satisfy energy conservation. 

In direct numerical simulation, we have calculated the potential energy U (e) [or, equivalently, 
the kinetic energy (de/dt) 2 ] as a function of e. As shown in FIG.O the decrease of U (e) agrees 



with our scaling and does not support the scaling U 



of Ref. & 
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Figure 6: Potential energy U{e) (where d e /L x = 0.01 and L y /L x = 47r in simulation) 



5 Discussions 

In this work, we have analytically elucidated the acceleration mechanism for collisionless re- 
connection driven by electron inertia. Let us interpret our result for tokamak parameters. For 
the m = 1 kink-tearing mode in tokamaks, r _1 = d e kB' corresponds to r _1 = d e q[ujAQ, 



where q[ is the derivative of the safety factor q at the q = 1 surface and ujao is the toroidal 
Alfven frequency at the magnetic axis. In order for the reconnection to be collisionless, the 
time scale r should be shorter than the electron-ion collision time r e = \i§d?Jr\, where 77 is the 
resistivity (at the q = 1 surface) and /x is magnetic permeability lfT3ll . For sample parameters, 
UM = 6.4 x 10V 1 , T e = 6keV, n = 3.5 x 10 19 m- 3 and q[ = 2.0m- 1 of TFTR (HQ, we obtain 
r = 90/xs and r e = 270/xs. Although the ratio r /r e can drastically change in proportion to 
r e _3 ^ 2 n 2 , these two time scales are not so separated but possibly similar in tokamak plasmas. 

Nevertheless, the time scale of explosion r predicted in this work is comparable to the ex- 
perimental sawtooth collapse times ~ 100/xs ifTSll . Note, inclusion of resistivity into Ohm's law 
® causes an additional decrease of the potential energy, one that would not prevent the release 
of free energy by inertia. In fact, our simulations exhibit nonlinear acceleration even with resis- 
tivity satisfying r /r e < 1. While the model used here is very simple, our result can be a central 
mechanism for sawtooth collapse. 

As might be expected, this explosive growth will be decelerated eventually before e reaches 
the equilibrium scale size L x (when the free energy starts to be exhausted). In tokamaks, we 
infer that the state of minimum potential energy is similar to the final state of Kadomtsev's 
model Q. But, if dissipation were sufficiently small, it would also corresponds to the state 
of maximum kinetic energy, where a strong convective flow remains. As shown in numerical 
simulations HH Q3), such a residual flow will cause a secondary reconnection and restore a 
magnetic field similar to the original equilibrium. 

We expect further applications of this variational approach to be fruitful for predicting strongly 
nonlinear and nonequilibrium dynamics of sawtooth collapses that other analytical methods fail 
to clarify. In addition to the theoretical estimation of the fast collapse time, a legitimate deriva- 
tion of a partial reconnection model (as well as associated loss of stored energy SW) would be 
made possible by extending the present analysis to more realistic two-fluid equations in toka- 
mak geometry. 
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